Sierra Nevada Time Series

Data preparation script for Sierra Nevada time series analysis

package ‘openxlsx’ was built under R version 3.6.2

Merge templates with 14C, C, and N data

Radiocarbon analyses for the 2001 samples were not run originally, but were completed on archived samples in 2020.

Fig. 16. Litter incubation \(\Delta\)¹⁴C-CO₂ (2019)

Caption: Mean \(\Delta\)¹⁴C-CO₂ for each site. Error bars show min and max of duplicate incubation samples. ^a) Data shown by site, without litter depth, ^b) Data shown by depth of litter layer, binned by climate zone.

2001 mean radiocarbon profiles

Fig. 16. Mean profile \(\Delta\)¹⁴C for 2001 samples

Caption: Mean \(\Delta\)¹⁴C by depth for each site in 2001. Error bars show ±1 standard deviation, solid vertical line shows \(\Delta\)¹⁴C of the atmosphere in the year of sampling.

2009 radiocarbon profiles

Fig. 16. Profile \(\Delta\)¹⁴C for 2009 samples

Caption: Profile \(\Delta\)¹⁴C by depth for each site in 2009. Solid vertical line shows \(\Delta\)¹⁴C of the atmosphere in the year of sampling. Error bars not shown as only a single replicate profile was analyzed per site.

2019 mean radiocarbon profiles

Fig. 16. Mean profile \(\Delta\)¹⁴C for 2019 samples

Caption: Mean \(\Delta\)¹⁴C by depth for each site in 2019. Error bars show ±1 standard deviation, solid vertical line shows \(\Delta\)¹⁴C of the atmosphere in the year of sampling.

Change in \(\Delta\)¹⁴C over time between 2001 and 2019

Fig. 16. Mean profile \(\Delta\)¹⁴C for 2001 and 2019 samples

Caption: Mean \(\Delta\)¹⁴C by depth for each site in 2001 and 2019. Error bars show ±1 standard deviation. Vertical lines show \(\Delta\)¹⁴C of the atmosphere in 2001 (solid) and 2019 (dashed).

Incubation \(\Delta\)¹⁴C-CO₂

Fig. 16. \(\Delta\)¹⁴C-CO₂ of 2019 bulk soil incubations

Caption: \(\Delta\)¹⁴CO₂ by depth for each site in 2019. One rep from GRrf 10-20 (the 10-20 cm increment sample from the cold granite site) is strongly depleted relative to the other rep: \(\Delta\)¹⁴C-CO₂ = -210.1. The highly depleted sample has been excluded for display reasons.

Fig. 16. \(\Delta\)¹⁴C-CO₂ of 2001 bulk soil incubations

Caption: \(\Delta\)¹⁴CO₂ by depth for each site in 2001. Note that some sites only have two depth increments. Similar to the 2019 dataset, one of the GRrf reps from the deepest depth increment was strongly depleted: \(\Delta\)¹⁴C-CO₂ = -469.1. Both points have been excluded for display reasons.

Fig. 16. \(\Delta\)¹⁴C-CO₂ of 2001 and 2019 bulk soil incubations

Caption: \(\Delta\)¹⁴CO₂ by depth for each site in 2001 and 2019. Different depth increments were sampled in 2001 and 2019. Points are the mean of laboratory duplicates; error bars are the measured values of each duplicate. Granite/cold point exlcuded for display reasons as it is strongly depleted.

Incubation vs. bulk soil \(\Delta\)¹⁴C

Fig. 16. Respiration data from incubations of 2019 and 2001 bulk soils.

Caption: Points show measured CO₂ production of laboratory duplicates as cumulative fluxes or daily flux rates by depth, lines show the means, and the ribbon represents the range.

Fig. 16. \(\Delta\)¹⁴C of 2019 bulk soil incubations and corresponding bulk soil

Caption: \(\Delta\)¹⁴C of bulk soil and respired CO₂ by depth for each site in 2019. Error bars show one standard deviation for bulk soil, points show mean of three replicate profiles for bulk soils and single observations for respired CO₂.

Fig. 16. \(\Delta\)¹⁴C of 2001 bulk soil incubations and corresponding bulk soil

Caption: \(\Delta\)¹⁴C of bulk soil and respired CO₂ by depth for each site in 2001. Points show mean of three replicate profiles for bulk soils and mean of laboratory duplicates for respired CO₂. The incubated soil samples are a composite made by homogenizing subsamples from each of the three replicate profile samples by depth. Error bars show one standard deviation for bulk soil and the measured values from laboratory duplicates of the incubated composite samples.

$emtrends
 PM d14c_mean.bulk.trend    SE df lower.CL upper.CL
 AN                0.514 0.144 42    0.224    0.804
 BS                0.873 0.169 42    0.532    1.214
 GR                0.952 0.136 42    0.679    1.226

Confidence level used: 0.95 

$contrasts
 contrast estimate    SE df t.ratio p.value
 AN - BS   -0.3593 0.222 42 -1.621  0.2481 
 AN - GR   -0.4387 0.198 42 -2.221  0.0792 
 BS - GR   -0.0794 0.216 42 -0.367  0.9287 

P value adjustment: tukey method for comparing a family of 3 estimates 

$emtrends
 eco d14c_mean.bulk.trend    SE df lower.CL upper.CL
 pp                 0.779 0.145 42    0.487    1.071
 rf                 0.777 0.317 42    0.137    1.417
 wf                 0.606 0.147 42    0.308    0.903

Confidence level used: 0.95 

$contrasts
 contrast estimate    SE df t.ratio p.value
 pp - rf   0.00222 0.349 42 0.006   1.0000 
 pp - wf   0.17373 0.206 42 0.841   0.6797 
 rf - wf   0.17150 0.350 42 0.491   0.8762 

P value adjustment: tukey method for comparing a family of 3 estimates 

Fig. 16. Regression of 2019 bulk soil incubations and corresponding bulk soil \(\Delta\)¹⁴C

Caption: Regressions of \(\Delta\)¹⁴C of bulk soil and respired CO₂ by depth for each site in 2019. Error bars show one standard deviation for bulk soil, points show mean of three replicate profiles for bulk soils and single observations for respired CO₂.

Time series: \(\Delta\)¹⁴C by depth (as measured)

Fig. 16. Time series of \(\Delta\)¹⁴C by depth, as measured

Caption: Points show mean of three profile replicates for 2001, 2009, and 2019 samples. Error bars show ± 1 standard deviation of the mean (only a single profile was analyzed in 2009). Stars show litter incubation \(\Delta\)¹⁴C-CO₂ for 2019 samples as a point of reference.

Spline fitting

Soils collected in both the 2001 and 2009 sampling campaigns were sampled by horizon, but the depth intervals differed between the two sampling years. In 2009, full profiles were excavated for each site, as opposed to the shorter profiles collected in 2001 from the GR and AN sites. Radiocarbon was measured on all three replicate profiles at each site for the 2001 samples, but only for one of the replicate profiles at each site in 2009, e.g. ANpp rep2, etc.

In order to compare the radiocarbon profiles between 2001, 2009, and 2019 we first interpolated both radiocarbon and carbon stock data at 1 cm intervals for each site in the datasets from each year. The carbon-stock-weighted radiocarbon values for any given target depth interval can then be calculated as a simple sum of the product of the carbon weight of each 1 cm increment (relative to the total carbon stock of the target depth interval) and its radiocarbon value. A monotonic cubic spline fit with Hyman filtering was used for the carbon stock interpolation (Wendt and Hauser 2013), and a mass-preserving spline was used to fit the radiocarbon data (Bishop, T.F.A., McBratney, A.B., Laslett, G.M., (1999) Modelling soil attribute depth functions with equal-area quadratic smoothing splines. Geoderma, 91(1-2): 27-45).

Error in `[.data.frame`(obj@horizons, , c(idcol, depthcols, var.name)) : 
  undefined columns selected

Fig. 16. Time series of bulk soil \(\Delta\)¹⁴C by 2001 depths (2001, 2009, 2019 samples)

Caption: Points for 2001 samples show the mean \(\Delta\)¹⁴C values at the measured depths. Points for 2009 and 2019 samples are spline-fitted estimates of \(\Delta\)¹⁴C predicted for the same depth intervals as measured in 2001. Error bars show ± 1 standard deviation of the mean of three replicate profiles for 2001 and 2019 samples (only a single profile was analyzed in 2009).

Fig. 16. Time series of bulk soil \(\Delta\)¹⁴C by depth (splined to 2019 depths)

Caption: Points for 2019 samples show the mean \(\Delta\)¹⁴C values at the measured depths. Points for 2001 and 2009 samples are spline-fitted estimates of \(\Delta\)¹⁴C predicted for the same depth intervals as measured in 2019. Error bars show ± 1 standard deviation of the mean of three replicate profiles for 2001 and 2019 samples (only a single profile was analyzed in 2009). NB: Only two depth intervals were measured at the cool and cold andesite sites (max depth of 27 and 28 cm, respectively), so linear extrapolation (using the slope of the last 1cm spline-fitted depth increment) was used to extend the profiles to 30 cm.

Fig. 16. Change in \(\Delta\)¹⁴C of bulk soil (panel a) and respired CO₂ (panel b) over time relative to the atmosphere

Caption: Points for 2019 samples show the mean \(\Delta\)¹⁴C values at the measured depths. Points for 2001 and 2009 (bulk only) samples are spline-fitted estimates of \(\Delta\)¹⁴C predicted for the same depth intervals as measured in 2019. Error bars for bulk samples in panel (a) show ± 1 standard deviation of the mean of three replicate profiles for 2001 and 2019 samples (only a single profile was analyzed in 2009); error bars for incubation samples in panel (b) show the values of the two reps, while the point represents the mean. NB: Only two depth intervals were measured at the cool and cold andesite sites (max depth of 27 and 28 cm, respectively), so linear extrapolation (using the slope of the last 1cm spline-fitted depth increment) was used to extend the profiles to 30 cm.

Initial modeling

The goal of this modeling exercise is to see how parent material and climate/ecosystem affect estimates of soil carbon ages and transit times. Bulk soil ¹⁴C observations from 2001, 2009, and 2019 will be used to constrain the carbon models, as well as observations of ¹⁴C-CO₂ from laboratory soil incubations of soils collected in 2001 and 2019. Previous work has indicated that the carbon stocks at these sites is likely at equilibrium, so we will apply the steady-state assumption to the modeling.

Two-pool models

One pool models have been shown repeatedly to be inadequate for describing soil carbon dynamics. However, as simple models are easier to constrain, we will start with a two-pool parallel and two-series models, as these are the simplest model system beyond the single pool approach.

The two-pool parallel model requires the following parameters: * decomposition constants for each pool (k₁, k₂) * input partitioning coefficient (\(\gamma\)) * steady-state carbon stocks (C) * inputs (I) * initial values of ¹⁴C 1 The two-pool series model requires the following parameters: * decomposition constants for each pool (k₁, k₂) * transfer coefficient (\(\alpha\)) * steady-state carbon stocks (C) * inputs (I) * initial values of ¹⁴C

Decomposition rates (k) are related to the amount of ¹⁴C in a pre-bomb system (fraction modern, F) at steady-state by the following equations (cf. Schuur, Druffle, and Trumbore, 2016): >Eq. 1

\[F = \frac{k}{k + \lambda}\] >Eq. 2

\[k = \frac{\lambda \cdot F}{1 - F}\] >where \(\lambda\) is the radioactive decay constant (1/8267).

As the decomposition rates will vary, the initial ¹⁴C content can be determined dynamically with equation 1.

Carbon stocks are known, while inputs will be estimated and are related to the steady-state conditions by the following equation: >Eq. 3

\[I = (k_{1} \cdot C_{1}) + (k_{2} \cdot C_{2})\] >where C₁ and C₂ are the carbon stocks of the two model pools.

Both stocks and inputs can be scaled to the known value of the total carbon pool once the steady-state parameters (k₁, k₂, and \(\gamma\) or \(\alpha\)) have been determined. Pool sizes are a function of the inputs and input partitioning coefficient at steady-state.

A Monte-Carlo Markov chain approach will be used for parameter estimation in combination with an initial optimization algorithm to determine the best set of initial parameters.

Workflow

Initial model fitting was performed for both model structures using generous parameter ranges [0, 1] for all three parameters (k₁, k₂, \(\gamma\) or \(\alpha\)). The initial parameter set was found by fitting the models by eye, followed by optimization with the function “modFit” (R package FME), using the Nelder-Mead algorithm. The best set of parameters found by modFit was then used as the input to a Monte Carlo Markov Chain (MCMC), using the function “modMCMC” (R package FME). The number of iterations for the MCMC optimization was set at 5000 intially, with delayed rejection employed to increase efficiency.

The sum of the mean squared error for the best parameter set was slightly lower for the parallel structure than for the series structure. Additionally, the overall mean error of the residuals was also lower for the parallel structure, moderately so for the bulk C observations but substantially so for the respiration observations (in andesite and granite soils in particular).

However, these initial fits yielded unrealistic parameter estimates for multiple sites, particularly at the lower depths. Additionally, the modFit output showed very high correlation between the parameters for both model structures (slightly higher for the two-pool series model).

Parameter optimization

Optimizing the parameter set requires imposing costs and optionally constraining the allowable range of values for each parameter. Given that we only have data for three time points, this is a relatively sparse data set for constraining these models. Accordingly, the optimization procedure will benefit from a priori constraints of the allowable parameter ranges. For example, since we assume that the system cannot be adequately modeled as a single homogenous reservoir, we will ensure that the optimization procedure cannot collapse the two-pool system into a single pool. This can be mitigated in the two-pool parallel optimization by constraining \(\gamma\) (i.e. the percentage of the inputs entering the fast pool) to a range of 50% to 95%. Similarly, for the two-pool series model structure we can constrain the range of the transfer coefficient to be between 0.0 and 0.1, ensuring that some carbon remains in the fast cycling pool.

Additionally, to enforce a relatively fast cycling pool and relatively slower cycling pool, we will loosely constrain the intrinsic decomposition rates as well (both model structures):

k₁: [0.02, 1.00] (50 to 1 year) k₂: [0.0001, 0.02] (10,000 to 50 years)

Finally, the models will be run to enforce steady-state, i.e. with unvarying carbon stocks. The amount of carbon observed in the system will be used in the cost function in addition to the radiocarbon observations made in 2001, 2009, and 2019. The inputs will be estimated from net ecosystem exchange (NEE) data measured at nearby eddy covariance sites: Blodgett experimental forest (AmeriFlux), Lower Teakettle (NEON), and Soaproot Saddle (NEON). Alternatively, using correlations between fluxes measured from these eddy covariance towers and GPP estimated from satellite retrievals of SIF, estimates can be made for inputs at the pixels corresponding to each site location.

Bayesian parameter estimation (MCMC)